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We study the free diffusion in two dimensions of active-Brownian swimmers subject to passive 
fluctuations on the translational motion and to active fluctuations on the rotational one. The 
Smoluchowski equation is derived from a Langevin-like model of active swimmers, and analytically 
solved in the long-time regime for arbitrary values of the Peclet number, this allows us to analyze 
the out-of-equilibrium evolution of the positions distribution of active particles at all time regimes. 
Explicit expressions for the mean-square displacement and for the kurtosis of the probability distri¬ 
bution function are presented, and the effects of persistence discussed. We show through Brownian 
dynamics simulations that our prescription for the mean-square displacement gives the exact time 
dependence at all times. The departure of the probability distribution from a Gaussian, measured 
by the kurtosis, is also analyzed both analytically and computationally. We find that for Peclet 
numbers < 0.1, the distance from Gaussian increases as ~ t~ 2 at short times, while it diminishes as 
~ t -1 in the asymptotic limit. 
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I. INTRODUCTION 

The study of self-propelled (active) particles moving 
at small scales has received much attention [1-9]. This 
is so since phenomena in nature like motion of plankton, 
viruses and bacteria, have an important role in many bi¬ 
ological processes as well as in industrial applications, 
and more generally, because active particles are suitable 
elements of analysis in the context on nonequilibrium sta¬ 
tistical physics. Moreover, confined matter made by self- 
propelled particles has been observed to behave in a very 
different way compared to matter conformed by passive 
ones, mainly due to the out-of-equilibrium nature of ac¬ 
tive systems. For example, it has been observed that 
active particles, confined in a box with sharp corners, ac¬ 
cumulate more at the corners of the container rather than 
spreading uniformly through the box. This effect origi¬ 
nates that the generated pressure on the walls do not be 
uniform along the container [10]. In contrast, matter con¬ 
formed by passive particles develop a uniform pressure in 
the container. It has also been reported that interacting 
active particles with different sizes and confined in a box, 
form clusters leaving empty regions inside the box [11], a 
very different scenario occurs with passive particles since 
these particles would spread homogeneously in the con¬ 
tainer. Thus we notice that the interplay among activity 
and confinement may provide novel properties to the lat¬ 
ter systems, hence the necessity of studying them. 

From a technological point of view, the study of ac¬ 
tive particles is also very relevant. To mention an ex¬ 
ample, the bioengineering community is constructing 
self-propelled micromachines [12 14] inspired by natu¬ 
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ral swimmers, and with the purpose of making devices 
able to deliver specialized drugs in precise region inside 
our body, or to serve as micromachines able to detect, 
and diagnose diseases [15-17]. These microrobots, in the 
same way as the smallest microorganisms, are subject 
to thermal fluctuations or Brownian motion, which is 
an important issue to take into account, since thermal 
fluctuations make the particles to lose their orientation, 
thus affecting the particles net displacement. Most clas¬ 
sical literature considering the effect of loss of orienta¬ 
tion due to Brownian motion and activity on the par¬ 
ticles diffusion, has used a Langevin approach, since it 
can be considered as a direct way of finding the parti¬ 
cles effective diffusion. Within this approach, isotropic 
self-propelled bodies subject to thermal forces [5, 18, 19] 
and anisotropic swimmers [18] have been treated in the 
absence and presence of external fields [20-23]. A Smolu¬ 
chowski approach to study the effective diffusion of ac¬ 
tive particles in slightly complex environments has also 
been undertaken. For example, Pedley [24] introduced 
a continuum model to calculate the probability density 
function for a diluted suspension of gyrotactic bacteria. 
Bearon and Pedley [25] studied chemotactic bacteria un¬ 
der shear flow and derived an advection-diffusion equa¬ 
tion for cell density. Enculescu and Stark [26] studied 
the sedimentation of active particles due to gravity and 
subject to translational and rotational diffusion. Simi¬ 
larly, Pototsky and Stark [27] followed a Smoluchowski 
approach and analyzed active Brownian particles in two- 
dimensional traps. Saintillan and Shelley [28] used a ki¬ 
netic theory to study pattern formation of suspensions of 
self-propelled particles. 

Self-propelled particles in nature move in general with 
a time-dependent swimming speed, like microorganisms 
that tend to relax (rest) for a moment and then to con¬ 
tinue swimming [29]. In this work, we assume that par¬ 
ticles move with a constant swimming speed, which is 
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a simplified model that has been supported by experi¬ 
mental work [20, 30] and has also been adopted when 
studying collective behavior [31]. Other more complex 
scenarios where for example active Brownian particles 
are subject to external flows, have been reported and 
solved following a Langevin approach [23]. 

In this paper we follow the approach of Smoluchowski 
and focus on the coarse-grained probability density of 
finding an active Brownian particle—that diffuses trans- 
lationally and rotationally in a two-dimensional, un¬ 
bounded space, and immersed in a steady fluid—at po¬ 
sition x at time t without actually knowing its direction 
of motion. We derive Smoluchowski’s equation for such 
probability density and we solve it analytically. 

We are able to use the Fourier analysis to translate the 
Fokker-Planck equation for active particles, into an infi¬ 
nite set of coupled ordinary differential equations with 
a hierarchy similar to the infinite BBGKY (Bogoliubov- 
Born-Green-Kirkwood-Yvon) hierarchy that appears in 
the Boltzmann theory, which in principle can be system¬ 
atically solved. We close the hierarchy at the second level 
and explicitly solve the first hierarchy equation, which 
corresponds to the Smoluchowski equation, and provides 
the probability density function (p.d.f) we are interested 
in. We transform back to the coordinate space, and we 
are able to explicitly find the p.d.f. of an active Brownian 
particle in the short and long time regimes. The mean 
square displacement (msd) for a swimmer at short and 
long times is also found, and classical results for the en¬ 
hanced diffusion of an active particle are recovered. We 
also find theoretical results for the kurtosis of the swim¬ 
mer, hence a discussion concerning the non-Gaussian be¬ 
havior at intermediate times of an active swimmer p.d.f. 
is also offered. Finally and with comparison purposes, 
Brownian dynamics simulations were also performed ob¬ 
taining an excellent agreement among our theoretical and 
computational results at short and long times for both, 
kurtosis and mean-square displacement of the swimmer. 


II. THE MODEL 

Consider a spherical particle of radius a, immersed 
in a fluid at fixed temperature T, that self-propels in 
a two-dimensional domain. The particle is subject to 
thermal fluctuations (modeled as white noise) in trans¬ 
lation and rotation, that is, £r(f) and £j?(t) respectively, 
where (£r> = (6?) = 0, {ki,T{t)£j, T (s)) = 2 D B S(t - s) 
and (£r(£)£r( s )) = 2Dn5(t — s). Here Db = k B T/dnrria 
and Dq are respectively, the translational and rotational 
diffusivities constants, with rj the viscosity of the fluid. 

The particle swimming velocity, JJ s (t), is written 
explicitly as U s (t)u(t), where we denote by u(t) = 
[cosy>(i),siny>(i)] (p(t) being the angle between the di¬ 
rection of motion and the horizontal axis) the instan¬ 
taneous unit vector in the direction of swimming, and 
U s (t) the instantaneous magnitude of the swimming ve¬ 
locity along u(t). Each of these quantities may be deter¬ 


mined from its own dynamics [32], and Langevin equa¬ 
tions for each may be written. Here we consider the case 
of a faster dynamics for the swimming speed such that 
U s (t) = Uq =const. In this way, the dynamics of this ac¬ 
tive Brownian particle, subject to passive-translational 
and active-rotational noise, is determined by its position 
x(t) and its direction of motion u(t), computed from 
p(t), that obey the following Langevin equations 


= U 0 u(t) + £ T (t), 

(la) 

= &(*)■ 

(lb) 


In this way, the temporal evolution of the particle’s po¬ 
sition [Eq.(la)] is thus determined by two independent 
stochastic effects, one that corresponds to translational 
fluctuations due to the environmental noise, the other, 
to the swimmer’s velocity whose orientation is subject 
to active fluctuations [Eq. (lb)]. These same equations 
have been considered in Ref. [33] to model the motion of 
spherical Platinum-silica Janus particles in a solution of 
water and ^ 02 - 

From now on, we use D ^ 1 and UqDq 1 as time and 
length scales respectively, such that Db = DbDq/Uq 
is the only free, dimensionless parameter of our analysis 
which coincides with the inverse of the so called Peclet 
number, which in this case corresponds to the ratio of 
the advective transport coefficient, Uq/Dq, to the trans¬ 
lational diffusion transport coefficient Db- For the Janus 
particles studied by Pallaci et al. in Ref. [34], we have 
that Dq 1 ~ 0.9 s -1 , Db ~ 0.34 /um 2 /s and swimming 
speeds between 0.3 and 3.3 fim/s. These values give Db 
between 0.035 and 4.2. Numerical results, on the other 
hand, are obtained by integrating Eqs. (1) using the 
Euler-Crome explicit scheme with a time-step of 0.005. 
In Fig. 1 seven single-particle trajectories are presented 
to show the effect of varying Db- The persistence effects 
are conspicuous for small values of Db (thick curves). 

An equation for the one particle probability density 
P(x , p, t) = (S(x—x(t))S(p—p(t))) can be derived by dif¬ 
ferentiating this P{x,ip,t) with respect to time, namely 

d 

-XT P{x, p, t) + U 0 u ■ VP{x , ip, t) = 

ot 

- - x(t))S(p - ip(t))} 

- V • (ir(t)S(x - x(t))S(ip - <p(t))), (2) 

where (•) denotes the average over translational and ro¬ 
tational noise realizations and V = (d/dx,d/dy). By us¬ 
ing Novikov’s theorem, we obtain the following Fokker- 
Planck equation for an active Brownian particle 

^-P{x, p,t) + U 0 u■ VP{x, ip, t) = 
ot 

02 

D B V 2 P{x,p,t) + Dn-^P{x,p,t). (3) 
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- £>„ = 0.001 

- D b = 0.00316 

— D a = 0.01 
- £>. = 0.0316 



FIG. 1. (Color online) Single particle trajectories for different 
values of D B = 0.001, 0.00316, 0.01, 0.0316, 0.1, 0.316, 1.0. 
Data is generated by solving Eq. (1) during 10 4 time steps, 
the initial position is chosen at the origin while the initial 
orientation is drawn from a uniform distribution in [0, 27 t). 

We now start to analytically solve Eq. (3). To do so, we 
apply the Fourier transform to Eq. (3) and we obtain 

^P(fc, p, t ) + iU 0 u ■ kP(k, p, t ) = 

d 2 ~ 

- D B k 2 P(k,p,t) + Dn-^-^P(k,<p,t), (4) 

where 


ficient of the expansion p n (k , t), namely 

^Pn = [e 2nDnt e~ ie p n -1 + 

e- 2nD ^e ie p n+1 \. (8) 

Note that we have introduced the quantities k x ± ik y = 
ke ±l9 . 2 

Our interest lies on Pq{x, t ) = / 0 ” dp P{x, p, t), which 
is given by the inverse Fourier transform of Po(k,t) = 
e -o B k as can b e checked from Eq. (7). Let 

us find po by solving Eq. (8) when the initial condition 
P n (k, 0) = p n (k, 0) = 8 nt o/27r is considered. This ini¬ 
tial condition corresponds to the case when the particle 
departs from the origin in a random direction of motion 
drawn from a uniform distribution in [0, 2n) and is equiv¬ 
alent to the initial condition P(x,p, 0) = (x)/2tt, 

where <5„ )m and S^^x) denote, respectively, the Kro- 
necker delta and the delta function. 

We first deal with the solution of Eq. (8) for long times. 
To do so, we consider only the first three Fourier coeffi¬ 
cients p n , i.e., those with n = 0, ±1, later on it will be 
shown that the other coefficients do not contribute in this 
regime. From Eq. (8) the three first Fourier coefficients 
satisfy 

^Po = -Y ik e~ Dnt [e~ l9 p- 1 + e l8 pi\ , (9a) 

j t P± i = -^ik [e^e^po + e -3 ™ e ±ie p ±2 ] , (9b) 

which after some algebra are combined into a single equa¬ 
tion for po, that is 


P(k, p, t) = (2n) 1 J d 2 x e lk x P(x,p,t), (5) 

and k = ( k x ,k y ) denotes the Fourier wave-vector. For 
Uq = 0, Eq. (4) has as solution the set of eigenfunc¬ 
tions {e - ( DBfe + D ^ n j w ith n an integer. Thus we 

expand P(k,p,t) on this set, expressly 


d 2 d„ C / 0 2 

dP m + Da dt m = -- kpn 


TJ 2 
u 0 
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- -^-k 2 e~ ADnt ( 


2i@~ | _—2 i6~ \ ( -i r\\ 

e p 2 +e P- 2 ) ■ (10) 


At this stage we argue that in the long time regime we 
may neglect the second term in the rhs of Eq. (10). This 
approximation leads to the telegrapher’s equation 


~ 1 

P(k,p,t) = — Pn(M)e- (DBfc2+% " 2) V^, (6) 

n=—oo 


with k = |fc|. This expansion allows to separate the ef¬ 
fects of translational diffusion contained in the prefactor 
e -D B k t f rom the effects of rotational diffusion due to 
active fluctuations. 

The coefficients of the expansion are given by 

p n (k,t) = e^ DBk2+D ^ t dpP(k,p,t)e~ in *. (7) 

Jo 

After substituting Eq. (6) in Eq. (4) and using th orthog¬ 
onality of the eigenfunctions we get the following set of 
coupled ordinary differential equations for the n-th coef¬ 


d 2 „ d_ Uq i2 ~ 

w ^ + D a - P „ = --k Pa , (11) 

which is rotationally symmetric in the fc-space and there¬ 
fore, giving rise to rotationally symmetric solutions in 
spatial coordinates if initial conditions with the same 
symmetry are chosen. One may check that Eq. (11) has 
the solution 


-sinojfet + costal , 

2 

( 12 ) 

withw^ = U'ok 2 /2 —Dq/4:. Therefore, by using this result 
one can explicitly write that Po(k,t) = e~ DBk t po(k,t), 
with e~ DBk * being the translational diffusion propa¬ 
gator. The solution in spatial coordinates, Po(x,t), is 


Po{k,t) = p 0 {k,Q)e Dnt/2 
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obtained after taking the inverse Fourier transform of 
P 0 (fc, t), which is given by the convolution of the transla¬ 
tional propagator with the probability distribution that 
retain the effects of active diffusion, i.e., 




P 0 (x,t) 


d 2 x'G (x — x', t)p 0 (x', t), (13) 


where po(x, t) is the inverse Fourier transform of po(k, t ) 
and 


G (x, t ) 


e ~x 2 /4D B t 

4irD B t 


(14) 


is the Gaussian propagator of translational diffusion in 
two dimensions. 

In the asymptotic limit (f —>• oo), at which the coher¬ 
ent wave-like behavior related to the second order time 
derivative in Eq. (11) can be neglected (mainly due to the 
random dispersion of the particles direction of motion), 
Po(x,t) tends to a Gaussian distribution with diffusion 
constant Uq/2Dq, namely 

e ~x 2 /4,(U 2 /2D n )t 

S- i,(UZ/2D a )C (15) 

Substitution of this last expression into Eq. (13) and per¬ 
forming the integral, one gets in the long-time regime 

e -x 2 /4(D B +U 2 /2Da)t 

F ° {x - t) = i^D7Twwm’ (16) 

from which the classical effective diffusion constant D = 
D B + Uq/2Dq is deduced [34]. 

In the short time regime Dot <C 1, we approximate 
Po(k,t) « po(k,0) = (2 tt)~ 1 and the p.d.f. in spatial 
coordinates corresponds to the Gaussian 


P 0 (x,t) 


l e ~x 2 /4D B t 

2n 4nD B t 


(17) 


Though it is a difficult task to obtain from Eq. (13) the 
explicit dependence on x and t of P 0 , a formula in terms 
of a series expansion of the operator V 2 applied to G{x, t) 
can be derived, to say 


FIG. 2. (Color online) Snapshots of the positions of 10 5 
particles taken at different times, namely from left to right, 
Dnt = 0.01, 0.1, 1.0, 10, and 100, for an inverse Peclet num¬ 
ber Db = 0.001. Notice the formation of a rotationally sym¬ 
metric ring-like structure (second panel from left) that devel¬ 
ops from a Gaussian distribution (far left panel) due to the 
effects of persistence in the motion of the particles. The struc¬ 
ture fades out as time passes (third and fourth panels from 
left) to reach a Gaussian distribution at long times (far right 
panel) . 


formation of a rotationally symmetric ring-like structure 
at the intermediate time Dot = 0.1 (second panel from 
left), mainly due to the persistence effects on the parti¬ 
cles’ motion. The structure fades out with time (third 
and fourth panels from left) and starts turning into the 
Gaussian distribution [Eq. (16)], as time passes (far right 
panel for Dot = 100). We must mention that in spite of 
the good agreement in the short and asymptotic limits 
given by (18), this solution does not provide an adequate 
description of the particles distribution in the interme¬ 
diate time regime for all values of D B , particularly for 
the small ones, for which the effects of persistence are 
conspicuously prominent. The reason of this failure lies, 
not in the general solution (13) but in the approximated 
equation (11), which corresponds to the two dimensional 
telegrapher’s. The solution to this equation [(12)] can 
not be interpreted as a p.d.f. since becomes negative at 
short times [35]. This undesired feature, results from the 
wake effect, characteristic of the solution of the two di¬ 
mensional wave equation [9]. Notwithstanding this and 
as is shown later on, the mean-square displacement com¬ 
puted from our approximation coincides with the exact 
one computed from numerical simulations, at all time 
regimes. 


III. MEAN-SQUARE DISPLACEMENT 


Po(x,t ) 


1 1 

1 P -Dnt/ 2U 1 

2tt ^ (2a)! 

s=0 v ' 



Pat 1 \ 

2 2s + 1 ) 


x 



G(x,t), 


(18) 


Note that the last formula can be rewritten in terms of 
d t instead of V 2 by using that V 2 G = D^dtG. 

In figure 2 snapshots of the positions of 10 5 particles 
taken at different times from numerical simulations, and 
for the inverse Peclet number D B = 0.001, are shown. 
The first panel from left {D^t = 0.01), corresponds to a 
p.d.f. close to the Gaussian given by Eq. (17). Note the 


What are the effects of self-propulsion on the diffusive 
behavior of the system is a question that has been ad¬ 
dressed in several experimental and theoretical studies, 
and the msd, being a measure of the covered space as 
a function of time by the random particle, has been of 
physical relevance in both contexts. Indeed, the msd is 
obtained in many experimental situations that consider 
active particles [18, 33, 36], hence a way of validating 
existing theoretical approaches. On the other hand, it is 
known that the diffusion coefficient of active particles de¬ 
pends on the particle density, due mainly to excluded vol¬ 
ume effects of the particles. Here we assume an enough 
diluted system to neglect those effects. 
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In what follows, an exact analytical expression for the 
msd will be obtained. From Eq. (11) one can show that 
in the coordinate space Po(x,t) satisfies 

S p ° +d 4 f »=(? +oA ) v ' fi 

+ D b ^ -D b V 2 ^ V 2 P 0 . (19) 

If last equation is multiplied by x 2 an integrated over the 
whole space, we obtain for the msd the equation 

^<® 2 (t)> + D Q j t (x 2 (t)) = 2 U 2 + 4D B D n , (20) 

whose solution can be easily found, namely 


{x 2 {t)) = 4^ (1-e Dnt ) 


2 (2 D B D n + U 2 ) 


n 2 

u o. 


[D n t-( l-e- Dnt )], (21) 


where we have used that ( d/dt ) (x 2 (t))\ t0 = 4 D B in 
this case. 

The linear dependence on time of the msd, expected 
in the Gaussian regime, is checked straightforwardly from 
Eq. (21). In the long-time regime we get 


(x 2 (t)) 


Dnt—too 


A 4 [D b + U ° 


2D n J 


( 22 ) 



FIG. 3. (Color online) Mean-square displacement in units of 
Uq /Dq as function of the dimensionless time Dnt for different 
values of the inverse Peclet number Db, namely 0.001, 0.0031, 
0.01, 0.0316, 0.1, 0.3162, 1.0. Solid lines correspond to the 
analytical expression given by Eq. (21), while the thin, dash- 
dotted lines correspond to the results obtained from numerical 
simulations (see text). 


We now use that 


(x 2 (t)) = -V 2 k P 0 (k,t) 


1 d d 
k dk dk 


P 0 (k,t) 


k—0 ’ 


(25) 

hence, at short times, the msd of the active particle is 
finally given by 


which is a classical result indicating that the effective 
diffusion of a self-propelled particle is enhanced by its 
activity [3]. In the opposite limit. 

(x 2 {t)) - >4 Dgt, (23) 

' ' Dnt-yO 

which once again the latter equation shows, the linear 
behavior with time of the msd. 


A. Comparison of the mean-square displacement 

In this subsection we compare the results obtained for 
the msd at short and long times using Eq. (19) with re¬ 
sults obtained following our new Fourier approach. 

For long times and with the explicit expression for 
P 0 (x,t) [Eq. (16)], we use the definition of the mean- 
square displacement, (x 2 (t)) = f d 2 x x 2 Pq(x, t) to easily 
recover Eq. (22). For short times, we expand Po(k,t) in 
powers of k, and we keep only the quadratic terms in k 
for the reason that will be clear immediately afterwards, 
we have 


Po{k,t) — 1 + 




t — 


^ P B Pnk 2 



( 24 ) 


(x 2 (t))^4D B t+(2D B D n + U 2 )t 2 , (26) 

which shows linear and quadratic behavior of the msd. In 
the limit t —>• 0, the above equation reduces to Eq. (23). 
Thus we have arrived to the same results by two distinct 
methods, namely, one method that obtained the explicit 
expression for the probability density and used it, and the 
second method that extracted information from Eq. (19) 
without solving it. 

In order to validate our theoretical findings, we have 
also performed Brownian dynamics simulations. Fig. 3 
shows a comparison among our theoretical prediction for 
any time given by Eq. (21) (solid lines), and our Brown¬ 
ian dynamics simulations (dashed lines) implemented to 
solve Eqs. (1) and averaging over 10 5 trajectories. Each 
plotted line corresponds to different values of D B . We ob¬ 
serve and excellent agreement among theory and Brow¬ 
nian simulations. Additionally, Fig. 3 indicates that a 
transition of the msd from a linear to a quadratic be¬ 
havior is not always the case. For large values oi D B 
the linear behavior over time dominates, whereas for D B 
smaller than 0.3, a transition starts to occur. The ab¬ 
sence of a transition from a linear to a quadratic behavior 
has also been reported by ten Hagen et al. [18]. They 
suggest that at short times, the msd of a self-propelled 
particle is affected by the initial orientation and the mag- 
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nitude of its force propulsion that make that a transition 
occur or not. 


where x T denotes the transpose of the vector x and E 
is the 2x2 matrix defined by the average of the dyadic 
product (x — (x)) T ■ (x — ( x )). In addition, it can be 
shown that for circularly symmetric distributions, as the 
ones considered in the present study, Eq. (27) reduces to 


IV. KURTOSIS 


K = 


> 4 (*))r 

(x 2 (t))r 


(28) 


Once we have obtained approximately Po(x,t) from 
Eq. (3), we wish to characterize its departure from a 
Gaussian p.d.f. as a function of time. This quantity 
has been measured experimentally is systems of Janus 
particles To this purpose we calculate its kurtosis k. given 
explicitly by [37] 

« = ([(*“ ( X )) T ^\ X - (®))] 2 ), (27) 

_I 


here (•),. stands for the radial average over the ra¬ 
dial probability density distribution, namely {’) ra d = 
/ 0 °° dr rP(r)(-). The analytical results obtained for the 
kurtosis will also be compared with those obtained from 
numerical simulations. 

As we can see, Eq. (28) requires the calculation of the 
fourth moment. After certain algebraic steps (see ap¬ 
pendix B), we explicitly find that the fourth moment is 
given by 


<* 4 (t)>r = Ae ~^ 2 


n 4 
U Q. 


-3 


Dnt (Dnt 


(l+4I> B (l + D B )) 

Dftt. ~ \ (Dfit 

3-^-(1 + 4 D b )+(-?- 


cosh I + 


(1 + 4D b (1 + D b )) 


**(¥ 


(29) 


The latter analytical expression is used in Fig. 4 to plot 
the kurtosis, to be precise, the ratio of Eq. (29) to the 
square of Eq. (21). In this figure, kurtosis is plotted as a 
function of the dimensionless time D&t for different val¬ 
ues of D b (dashed lines), while the corresponding exact 
results from Brownian dynamics simulations are shown 
in symbols. Note that for these last ones, 4 < n < 8. 

We observe that the Jcurtosis shows a clear non¬ 
monotonic behavior for D B < 1, namely, it starts to 
diminish from k = 8 as time passes, until it reaches a 
minimum at a time around D^ 1 that depends on the in¬ 
verse of the Peclet number D B . This behavior has been 
observed in experiments with Janus particles in three di¬ 
mensions [33] and theoretically predicted in one quasi- 
one-dimensional channels [38]. In the limit of vanishing 
translational diffusion, the kurtosis is a monotonic in¬ 
creasing function of time, with n = 4 as its minimum 
value [9]. At later times the kurtosis starts to increase 
reaching the value 8 in the asymptotic limit. This be¬ 
havior becomes less and less evident as the translational 
diffusion coefficient surpasses the effective diffusion coef¬ 
ficient that originates in the rotational diffusion, Uq/Dq. 
Due to this effect, we note a better agreement between 
our analytical approximation and our simulation results, 
as D b is increased. On the other hand, the persistence 
effects induced by orientational correlations are revealed 
in the diminishing regime of the kurtosis, as is shown in 
Fig. 4 for the cases in which D B < 1 (see Table I for 


the values of n corresponding to the snapshots shown in 
Fig. 2). In this regime the orientational correlations sur¬ 
pass the effects of translational noise hence allowing the 
particles to move persistently outwardly from the origin 
and giving rise to circularly-symmetric, ring-like distri¬ 
butions as the one shown in Fig. 2 (second from left) 
at Dnt = 0.1 for D B = 0.001. Afterwards, this ring¬ 
like distributions turns into a Gaussian distribution as 
the orientational correlations die out. Notice that the 
discrepancy between the results given by our analytical 
approximation and the exact results lies on the wave-like 
effects given by the second-order-time-derivative term of 
Eq. (11) which, as shown in [9], gives 8/3 ~ 2.6667 as 
the value for the kurtosis in the vanishing translational 
diffusion. For D B = 0.001, the minimum of the exact 
kurtosis is about 4.1352 at t « 0.39 -Dq 1 while from our 
analytical approximate n « 3.0184 at 0.239 -Dq 1 . 

In the next subsections we find asymptotic expression 
for the kurtosis for both short and long times that confirm 
the Gaussianity of the distribution function for these time 
regimes. 


A. Kurtosis at short times 

In ordet to find k for t 0, we expand Po(k,t) in a 
power series in k, use the fact that = U^k 2 /2 — Dq/ 4, 
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TABLE I. Kurtosis valuesat times Dnt that are close to those 
values chosen in Fig. 2 ( D b = 0.001). 


Dnt 

0.01 

0.1 

1.01 

10.245 

100.08 

K 

5.9598 

4.3111 

4.2729 

6.7621 

7.8713 



B. Kurtosis at long times 


In ordet to find k for t —>■ oo, we proceed as before, 
we expand PP 0 (k,t) in power series of k. use that = 
Uq k 2 /2 — Dq/ 4, but now keep all the terms in the series 
expansion. After applying Eq. (Bl) we get that in the 
long time limit 




2 5 U 4 r, , 4 D b Dq 

D 4 1 + U 2 

u a 1 




2 



2 


(32) 


Finally we use Eq. (22) to get for long times that 


k = 


(x 


Kt)) 


4 ( D b + 


_GL) 

2 Dq, J 



(33) 


which once again, Eq. (33) shows that kurtosis has a 
Gaussian behavior for long times. This result has also 
been validated using Brownian dynamics and shown in 
Fig. 4, where for long values of the dimensionless time, k 
tends again to a value of 8. 


FIG. 4. (Color online) Kurtosis of the particles distribution 
as function of the dimensionless time Dnt for different val¬ 
ues of Db, namely 0.001, 0.0031, 0.01, 0.0316, 0.1, 0.3162, 
1.0. Dashed lines correspond to the exact analytical expres¬ 
sion given by the quotient of Eq. (29) and Eq. (21), while 
the symbols mark the values of the kurtosis obtained from 
numerical simulations. Thin-dotted lines mark the values 8, 
4, 8/3, that correspond respectively in two dimensions to: a 
Gaussian distribution, a ring-like distribution and wave-like 
distribution [9]. 


and keep terms only of order k 4 to get 


P 0 (k,t) ss k z 


+ 


Dq C 4 t 5 C 4 t 4 Dq C 2 D B t 4 

~2 5!~~ + “iT + ~2 3! 

c 2 D B t 3 , DqD b D D 2 B t 2 


+ 


= k 4 f(t). 


(30) 


We now use Eq. (Bl), that together with Eq. (26), one 
finally gets for short times that 


V. THE ACTIVE FLUX: THE VELOCITY 
FIELD V(a:,f) OF ACTIVE PARTICLES 

The Fourier expansion given by Eq. (6) allows us to 
obtain the hydrodynamic fields. At the moment we 
have explicitly calculated Po(k,t) that corresponds to 
the Fourier transform of the density field p(x,t). The 
Fourier transform of the components of the velocity field 
V{x,t) are related to p±i(k,t) in the following way. By 
definition 

dip cos ipP{k, ip, t), (34a) 

dp siiupP{k,ip, t), (34b) 

and after substitution of Eq. (6) in Eq. (34) we get 

V x (k,t) =-y e -(D B k 2 +D n )t ^k,t) +pl(-k,t )\, 

(35a) 

V y (k,t) =y e ~ (DBk2+Dn)t K(-M) ~Pi{k,t )\, 

(35b) 


V x (k,t) = U 0 


V y (k,t) = U 0 


K 


_ 4 3 f(t) _ = g 

(4 D B t + 2 D b DqD + U 2 t 2 f 


(31) 


which means that our derived p.d.f. has a Gaussian be¬ 
havior at short times. This result has been validated 
numerically and illustrated in Fig. 4, where kurtosis is 
plotted versus time. It clearly can be seen that for short 
times k has a value of 8. 


where we have used that p- n (k,t) = p^(—k,t). With 
these expressions, it is straightforward to verify that Eq. 
(9a) corresponds to the continuity equation dtp + V • 
V(a;, t) = 0. At the long time regime, where only the first 
three Fourier modes are needed, we have that p\{k, t) can 
be computed from Eq. (12) if we neglect the term that 
involves P 2 , i.e., from 

±Pi = -^i ke -" e D ‘' t Po. 


(36) 































We now use the explicit expression for p 0 , that is Eq. (12), 
with po(k, 0) = 1/(27t) to integrate Eq. (36). After some 
steps we find that 

Pi(k,t) = ( ik x + k v ) e Dut / 2 Sm0Jkt . (37) 

47T UJk 

If we insert the latter into Eqs. (35a)-(35b) we obtain 

vjk.t) = 

47T Wfe 

Vy(k,t) = 

47T UJk 

After inverting the Fourier transform, this pair of expres¬ 
sions can be brought into the form 

y(x,t) = -|vF( a ;,t), (38) 

where the function F(x,t) depends on To, explicitly 

F(x,t) = J dse Da( - t ~ s ^ J d 2 x'G(x — x', t — s)Pq(x', s). 

( 39 ) 

Last expression corresponds to a non-Fickian constitu¬ 
tive relation that considers nonlocal effect in space and 
time. This relation, together with the continuity equa¬ 
tion, gives rise to a non-local (in both, space and time) 
diffusion-like equation. 

In a similar way, one can systematically find higher 
Fourier modes for the probability density function, 
P(fc, ip, t), in the wave vector space. 

VI. FINAL COMMENTS AND CONCLUSIONS 

Though Eq. (21) has been derived from the long time 
approximation, it does give the correct time dependence 
at all times as it is shown in Fig. 3, where our analyt¬ 
ical formula (solid lines) is compared against numerical 
simulations (dashed lines). This fact can be understood 
directly from Eq. (10), since it can be shown that the 
inhomogeneous term does not contribute to the mean 
squared displacement (see appendix A), therefore agree¬ 
ing with the expression obtained from the telegrapher’s 
Eq. (11). 

We want to point out that though the linear time de¬ 
pendence of the msd, characteristic of normal diffusion, is 
reached at times D^t ~ 1 (see Fig. 3), the Gaussian be¬ 
havior of the distribution is reached at much later times 
(see Fig. 4), suggesting that at finite times, D^t > 1, 
normal diffusion does strictly corresponds to Gaussian 
distributions. The kurtosis computed from our approxi¬ 
mation shows that it tends asymptotically to Gaussian, 
k = 8, as (8 — Z?nt) _1 . In contrast, in the short time 
regime, the distribution departs from a Gaussian behav¬ 
ior with time as (8 — Do,t)~ 2 (see Fig. 5). 



FIG. 5. (Color online) The departure from Gaussian behavior 
of Po (x,t), measured by 8 — k, differs in the short-time regime 
from the long-time regime. In the first case the departure 
captured by the power law (Dnt)~ 2 while in the long-time 
regime such departure is well described simply by (T>nt) _1 . 


In summary, we have found a method and obtained 
with it, an analytical solution to the Smoluchowski equa¬ 
tion of an active Brownian particle self-propelling with a 
constant velocity. We used a Fourier approach and ex¬ 
ploited the circular symmetry of the distribution func¬ 
tion, to obtain an infinite system of coupled ordinary dif¬ 
ferential equations for the Fourier modes of the probabil¬ 
ity density. Our formalism showed that in order to have a 
whole description for the particles diffusion, we only need 
the first Fourier coefficient, since it was found that higher 
Fourier modes do not contribute to the mean-square dis¬ 
placement of the particle. With the explicit p.d.f. in 
hand, we calculated the particle mean-square displace¬ 
ment for both long and short times, recovering classical 
results of enhanced diffusion due to activity. We also val¬ 
idated these findings by performing Brownian dynamics 
simulations that showed and excellent agreement among 
theory and simulations. In addition, we also performed 
an asymptotic analysis for short and long times, to the 
forth moment of our p.d.f. This calculation enabled us 
to have analytical results for the kurtosis. The asymp¬ 
totic results for this measure were validated with Brow¬ 
nian dynamics simulations that showed that the p.d.f. of 
an active particle always starts with a Gaussian behav¬ 
ior, and depending on the vale of the ratio Db, it may 
evolve to a ring-like distribution, but it always returns 
to a Gaussian distribution for long times. Finally this 
work shed light on mathematical procedures to solve a 
Smoluchowski equation for an active Brownian particle. 

Future venues for this research, would be to solve an¬ 
alytically a Smoluchowski equation for confined active 
particles and/or particles interacting among them. 
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and therefore, using Eqs. (A1) we finally have that 

+ Dnj t (x 2 (t)} 0 = 2 U 2 + AD b Dq, (A4) 

which coincides with Eq. (20), deduced from the tele¬ 
grapher Eq. (11) when considering only the first three 
Fourier modes po ,P±i- 


Appendix A: Mean Square Displacement 


Appendix B: Calculation of the Fourth moment of 

Po{x, t) 


If we apply the operator — V| to the relation Po(k, t) = 
e -D B k tj and evaluate at k = 0, one gets that 

{x 2 (t)) = AD B t + (x 2 (t))o, which leads straightforwardly 
to the pair of relations 


The fourth moment ( x 4 (t)) r can be obtained through 
the prescription 


(x 4 {t)) 


lAfcAV 

k dk dk ) 


Po{k,t) 


(Bl) 


^(® 2 (t)) = 4 D B + ^(x 2 (t)) c 


(Ala) 

(Alb) 


Though it is possible to use this prescription directly on 
Po(k,t), we prefer to apply it to the series expansion in 
powers of k of Po(k,t), this is due to the complicated 
dependence of Pq on k. Such expansion is given by 


where (-)o denotes the average of (•) with p 0 (x,t) as the 
probability measure. The quantity -^(x 2 (t ))o can be 
computed directly by applying — V| to Eq. (10) and eval¬ 
uating at k = 0, which results in 

-j^(x 2 (t))o+Dn — (x 2 (t)) 0 = [V):(fe(, + k 2 )p 0 \ fc=Q 

+ ^-e iDnt { V^, [( k x — ik y ) 2 p-2 + ( k x + ik y ) 2 p2\ } fc=0 

(A2) 

It is easy to show that the first term of rhs of the last 
equation is 2Uq while the next term, that involves the 
modes p ±2 vanishes identically, thus we get 

^(x 2 (t)) 0 + D n ^(x 2 (t))o=2U$, (A3) 


P 0 {k,t) = e 


_ p ~Dat/2 


00 m 

E E 

n,m=0 j =0 


m\ {—D B t) n ( it) 2m 
j) n! (2m)! 


x 



Dnt \ 

2(2m + 1) J 



(B2) 


where we have used the explicit dependence on /c, of Uk, 
given in section II. 

Since 


(is4) 2 ‘ 2l ” +i> =[ 2 <" + « i 2x 

[2(n + j — l)] 2 k 2 ( n+ i~ 2 \ (B3) 

the only terms that survives when evaluating at k = 0 
are those for which n + j =2 with the condition that 
0 < j <m. Thus we can write 




2m+2 


m —0 


(¥)' 


1 + 


Dnt 1 \ (m. + 2 )(to + 1) 


2 2?n + 5 / 2(2?n + 4)! 




+ 


Dnt 1 


m + 1 


2 2 m + 3/ (2m. + 2)! 


■Dill 


Dnt 1 
2 2m+lj2(2m)\ 


(B4) 


The sums can be written in terms of hyperbolic function, and after some algebraic steps Eq. (29) is recovered. 
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